NJASA- CHr n^QQ, 


NASA Contractor Report 178126 
ICASE REPORT NO. 86-36 

ICASE 


NASA-CR-178126 

19860022717 


A 


STABILITY ANALYSIS OF INTERMEDIATE BOUNDARY 
CONDITIONS IN APPROXIMATE FACTORIZATION SCHEMES 


Jerry C. South, Jr. 
Mohammed M. Hafez 
David Gottlieb 


Contract Nos. NAS1-17070, NAS1-18107 
May 1986 


INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING 
NASA Langley Research Center, Hampton, Virginia 23665 


Operated by the Universities Space Research Association 


NASA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton. Virginia 23665 


LIBRARY COPY 

St! 3 8 i986 



LANGLEY RESEARCH CENTER 
LIBRARY, NASA 
HAMPTON, VIRGINIA 



STABILITY ANALYSIS OF INTERMEDIATE BOUNDARY CONDITIONS 


IN APPROXIMATE FACTORIZATION SCHEMES 


Jerry C. South, Jr. 

NASA Langley Research Center 


Mohamed M. Hafez 
University of California, Davis 


David Gottlieb 
Brown University 


Dedicated to Milton E. Rose 
on Occasion of his 60th Birthday 


Abstract 

The paper discusses the role of the intermediate boundary condition in 
the AF2 scheme used by Holst for simulation of the transonic full potential 
equation. We show that the treatment suggested by Holst led to a restriction 
on the time step and suggest ways to overcome this restriction. The 
discussion is based on the theory developed by Gustafsson, Kreiss, and 
Sundstrom and also on the von Neumann method. 
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INTRODUCTION 


Approximate factorization schemes are widely used to obtain efficient 
solutions to problems in Computational Fluid Dynamics. In many cases, 
they have provided a significant Increase in efficiency over previously-used 
solution methods in particular problems. Some outstanding examples are the 
classical Alternating-Direction- Implicit method of Peaceman and Rachford [1], 
the Briley-McDonald Linearized Block Implicit scheme [2], and the Beam and 
Warming [3] Approximate Factorization (AF) scheme for the compressible Navier- 
Stokes equations. In the transonic potential-flow area, some AF schemes which 
have significantly improved solution efficiency are the work of Ballhaus and 
Steger [4], Ballhaus et al. [5], Holst [6], [7], and Jameson [8]. 

All of these schemes have the common feature that the solution procedure 
is broken down into a sequence of easily-implemented stages; i.e., easily- 
inverted matrix factors. Each of the stages usually requires boundary 
conditions for an "intermediate" variable (vector) which is not always a 
consistent approximation to the solution function desired. This feature can 
make satisfaction of implicit boundary conditions difficult, at best, and 
impossible, at worst. Dwoyer and Thames [9] demonstrated serious boundary- 
condition problems associated with the class of AF schemes called "Locally 
One-Dimensional," even in explicit schemes. 

The present paper further highlights the importance of intermediate 
boundary conditions by focusing on a specific example — a boundary-induced 
stability restriction in Holst's AF2 scheme [6] for the transonic full- 
potential equation. An analysis of the effect of the intermediate boundary 
condition is given by use of the usual von Neumann method and also the methods 
of Gustaffson, Kreiss, and Sundstrom [10] and Osher [11]. 
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ANALYSIS 


Holst's scheme is a variation of the AF2 schemes presented in References 
4 and 5. It will be referred to herein as "AF2Y," since in its implementation 
the y-operator is split, rather than splitting the x-operator as in References 
4 and 5. For the purpose of analyzing the intermediate boundary-condition 
problem, it is illuminating to study the application of AF2Y to the two- 
dimensional (2-D) Laplace's equation in a rectangle. The present analysis is 
valid only for the subsonic flow condition, which is simpler by far than the 
transonic case. However, it is reasonable to assume that if boundary-induced 
instability is present in the subsonic case, it will also occur in the 
transonic case. In practice this was true. 


The Discrete Problem 

The following thin-airfoil problem is thus considered: We wish to solve 

the Laplace difference equation for the disturbance velocity potential 



(a6 + 

xx 


b6 )$.. = 0 

yy jk 


(i) 


where a and b are constant coefficients and 
difference operators; e.g.. 


6 

xx 


and 


6 are central 

yy 


6 xx <J> jk ^j+l.k 



+ # 




( 2 ) 


The boundary conditions are set on a rectangular region with Dirichlet 
conditions, <p = 0, set on three sides (left, top, and right), representing 
vanishing disturbances, and a Neumann condition at the bottom boundary, 
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representing a thin-airfoil f low-tangency condition: 


4> y = s(x) 


at y = 0. 


A discrete analog of Eq. (3) at k = 1 can be written as: 


($ + $)*.,= 2Ays(x) 

y y 


where we use the following notation for one-sided, two-point differences: 


( 3 ) 


(4) 


Vjk ~ *j,krU " *jk 

Vjk = *jk ' +j,k-l- 

The difference operator (1) requires evaluation of 
k = 1. Since this operator can be written as: 


(5) 

( 6 ) 

6 <f> . . at the boundary 

yy J , 1 


6 


yy 



(7) 


Equation (4) is used to eliminate t <|k i» which calls for a value of 
below the boundary k = 1 . Thus, the difference operator at k = 1 is: 



L B <t>j,i = + 2b ^ y ^j,l " 2b 4ys( x )- 


( 8 ) 
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The AF2Y Scheme 

The AF2Y scheme models a hyperbolic equation, 
an iteration scheme: 


o<f> yt = V 2 <j>, 


and is used as 


(a + b,J )(- ab„$ - a6 )Ad>., = auLij) 1 ? 

1 y 2 y xx T jk T jk 

where n is the iteration counter, 


(9) 


b l b 2 " 


( 10 ) 


and A<|> is the correction 


■ *£' - *£• < u) 

The scheme is implemented in two stages: 

(a + b 1 J y )f jk = ati)L(j>j k (12) 

(- ab 0 £ - a6 )A<b.. = f .. . (13) 

2 y xx Y jk jk 

The intermediate variable f is defined by Eq. (13). The parameter a 
corresponds to a reciprocal "time" step, At *, and is usually cycled between 
small and large values to obtain rapid convergence. The parameter w 
corresponds roughly to a relaxation factor which is usually close to 2. 

The first stage (12) is bidiagonal, proceeding from the bottom boundary, 
k = 1, to the last interior row of mesh points, k = K - 1, for every j. The 
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second stage (13) is a tridiagonal solution which proceeds row-by-row, from 

k=K-l to k = 1, to obtain the correction The second stage is 

initiated with the condition A4> . v = 0, corresponding to the vanishing 

J ,k. 

disturbance, 4> = 0, at k = K. 

The Intermediate Boundary Condition 

The main problem in implementing the scheme is how to initiate the 
bidiagonal solution for f at k = 1. It seems reasonable, at first sight, 
to use a derivative condition on f at the boundary, as Holst [6] did; i.e., 

t f. . = 0. (14) 

y i»i 

Comparison of Eqs. (14) and (12) implies that 





(15) 


If this procedure is used with no further modification, it is unstable for 
small values of a (or large "time" steps) and fixed w as described next. 

Stability Analysis 

A von Neumann (VN) analysis shows that the interior scheme (9) is stable 
for all modes under the restrictions 

0 < to < 2 (16) 


a > 0 


(17) 
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However, the boundary scheme, implied by Eqs. (15) and (13) taken together, is 
another matter. 

A boundary condition more general than Eq. (14) for f can be 
considered. Let a "dummy-point" value for f be given as: 


f . n = yf . . . 
J,0 J,1 


(18) 


Then the equation for fj^ is, from Eq. (12), 


(a + bl )f. }1 = cu-L^j + Tbjf^j 


(19) 


and Eq. (13) yields: 


f j,l ' “ b 2 ? y ' * S « )A b,l = 


( 20 ) 


where 


n = 


am 


a+bj (1-y) 


( 21 ) 


To carry out a VN analysis, we substitute into Eq. (20) trial solutions 


i n = G n e i(jpAx+kqAy) 


( 22 ) 


where i = / -1 , p and q are wave numbers, and G is the amplification 

factor, to obtain: 


(aB + 2Ab x - iaE) (G -!)=-' 2nb 1 (A + B - iE) 


(23) 



where 


A = a (1 - cos ?) 
B = b (1 - cos n) 
E = b sin n 
£ = pAx 
n = qAy 



(24) 


The stability condition, 


< 1 , 


reduces to: 


n{(A + B)[ (2 - nJbj^A + (a - ftb^B] + (a - Gb^E 2 } > 0. (25) 


To maintain the inequality (25) the following stability restrictions are 
easily deduced: 

0 < fl < 2 (26) 

a > b l fi. (27) 

For the case y = 1 , corresponding to the backward-Neumann condition 
on f (Eq. (14)), restrictions (26) and (27) reduce to Eq. (16) and 


a > b^ to (y = 1). 


(28) 


The restriction (27) enforces a "time" step limitation on the scheme for fixed 
ft, which will slow convergence; or a reduction in ft, according to: 
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which in fact yields fast convergence and ensures stability. 

It is noted that another useful type of boundary condition for 
given by 



-Mr 


n 

3 ,1 


f is 


(30) 


which gives the same form as Eq. (20) for fj ^ with 


a(u) + 8) , „ . . 

° - a + b° ' (31 > 

Both classes of schemes are implemented by initiating the bidiagonal march for 
f using Eq. (20), under restriction (29). 

The restriction (29) was verified numerically in both a constant- 
coefficient, Cartesian-coordinate computer code for Laplace's equation and in 
the "TAIR" code [12] by using fixed values for a (i.e., no a-cycling) and 0, 
and for various values of the coefficient bj. In all cases, convergence was 
obtained when the restriction (29) was obeyed; and divergence occurred when it 
was violated. 

The experiments with the TAIR code were especially interesting, since the 
coefficient bj varies along the airfoil surface. The test case chosen was 
the default "0"-type mesh for an NACA 0012 airfoil. It was found that the 
arithmetic mean of bj along the surface presented the crucial condition, 
rather than the maximum value, as might be expected. 

The question arises as to why the TAIR code, which implements the AF2Y 
scheme with the boundary condition (14), operates so well since a is cycled 
between small values, which violate the restriction (28) and large values. 
The answer seems to be that a is increased within several rows adjacent to 
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the boundary to a value which (in the default mesh) meets the restriction 
(28), when smaller values of a are used in the remaining interior field. 
This "fix" was developed empirically by the authors of Reference 12; without 
this fix the code diverges. This procedure is not recommended in general, 
since it requires a discontinuous change in a. The assumption in the 
development of the factored scheme (9) is that a is constant throughout the 
mesh. 

A seemingly attractive scheme, involving a discontinuity in a at the 
boundary, is as follows: Initiate the solution for f using Eq. (15) with 

a) = 1 , and change the second stage (13) at the boundary to: 


f-2b£ - a6 ) A*. , = f . . = <j>. . , 

y xx ; *3,1 j,l B *3,1 


(32) 


This procedure exactly annihilates the boundary residual (in the linear case) 
and represents a fully implicit satisfaction of the surface boundary 

condition. However, the factored operator at line k = 2 is no longer the 

interior-point operator, since the term -abj^ in the inner factor is 

changed to -2b£^ discontinuously. It is possible to analyze such a scheme 

by the methods presented herein, but the line k = 2 must be considered as 

part of the boundary scheme. No details will be given here, but the analysis 
shows that setting u> < 4/3 at k = 2 guarantees linear stability of the 
overall scheme. However, the amplification factor modulus | G | exceeds unity 
only in a narrow frequency range of small n (Eq. (24)) when w > 4/3. 

Numerical experiments; showed no sensitivity to the value of to at k = 2. 
This scheme was always stable in tests with a constant-coefficient Cartesian- 

mesh code, even with m = 1.8 at k = 2. If the scheme was unstable for 
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highly stretched grids, setting w < 4/3 at k = 2 did not stabilize 

the scheme. In the variable-coefficient, nonlinear case, such a scheme is no 
faster than, and not as robust as, the scheme (20) with restriction (29). 

Review of the Stability Theory 

It is well known that in general the von Neumann analysis at a single 
line is neither sufficient nor necessary for checking stability. Trapp and 
Ramshaw [13] pointed out the usefulness of the VN analysis to study boundary 
schemes but recognized that no theoretical justification was known. 

We wish to review briefly the stability theory for finite-difference 
approximations to initial boundary-value problems. A necessary condition for 
the stability of such a scheme is the Ryabenkii-Godunov condition. It states 
that the numerical scheme is unstable if there exists a solution of the type 

^ >k = G n l G l > 1 (33) 

for the inner scheme and the boundary scheme. (It is also sufficient to check 

one boundary at a time.) Substituting (33) into (12) and (13) one finds that 

ip . , satisfies a constant coefficient second-order difference scheme whose 
J t k 

solution is 

* . . - U k e 1(jp6x) . (34) 

3 > K 

Actually there are two possible p's, but it is readily verified that only one 
of them satisfies |p| <1 for |g| > 1; and, therefore, it is not a valid 
solution for the quarter-plane problem. 
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In Appendix A we show that there exists a solution of the form (34) to 
(12), (13), and (20) such that j G j > 1 and |p| < 1 if (29) is not 

satisfied. This proves that the scheme is unstable. By instability here we 
mean that unbounded solution occurs after a fixed number of time steps for any 
mesh — it precludes the possibility of reaching steady state. 

It should be noted here that VN analysis of the boundary scheme does not 
predict the existence of solutions of the form (33) with j G | > 1. In fact, 
Gottlieb and Turkel [15] gave an example of a boundary scheme (Scheme VI, p. 
184 of Reference 15) coupled with a variant of MacCormack's scheme in the 
interior which is conditionally stable, yet the VN analysis of the boundary 
scheme shows unconditional instability. However, Goldberg and Tadmor showed 
that for a dissipative interior scheme (i.e., amplification — factor modulus 
bounded away from unity for all nonzero modes) VN stability of the boundary 
scheme excludes the possibility of an eigenvalue or a generalized eigenvalue. 
By an eigenvalue we mean a solution of the form (34) with |g| > 1 whereas a 
generalized eigenvalue is G such that |g| = 1. Thus, if the condition 
stated in (29) is satisfied no eigenvalue or generalized eigenvalue exists. 
In Appendix A we show it directly. The theory of Gustafsson, Kreiss, and 
Sundstrom [10] (see also, Osher [11]) states that for a system of first-order 
hyperbolic equations stability is assured if there is no eigenvalue or 
generalized eigenvalue. While their theory does not apply directly to the 
equation 

0 Py t = V 2 P, 


it can be modified to include this case. 
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As a concluding remark we should note that stability here implies 
convergence in the sense of Lax — the numerical solution converges to the 
analytic one as the mesh size tends to zero for fixed time t. This is 
clearly only a necessary requirement to reach steady state. 

Two-Dimensional Numerical Results 

A limited number of numerical tests for cases involving stretched grids 
and nonlinear transonic flow have convinced us that the discontinuous-a 
schemes (e.g., Eq. (32)) are not as reliable as the scheme using Eq. (20) with 
restriction (29). Some numerical results are presented in Tables 1 and 2. In 
the tables, the following identification is used for the various boundary 
schemes : 

Scheme I: Original TAIR scheme; Eqs. (13) and (15) at boundary, with a 

increased at 3 lines adjacent to boundary to satisfy restriction 
(28) with 10% safety margin. 

Scheme II: Exact annihilation of boundary residual; Eqs. (15) and (32), 

with to = 1 at boundary only. 

Scheme III: Eq. (20) and restriction (29) with 10% safety margin. 

Table 1 shows a series of numerical tests for incompressible flow over a 
circle, with varying degrees of mesh stretching near the boundary. The TAIR 
code was used with u = 1.8 at all points except as noted in schemes II and 
III, and with the default settings for the a-cycle (a m ^ n = 0.07, a max = 
1.5). The mesh contained 101 points uniformly spaced around the circle and 21 
points in the radial direction with stretched spacing. The first column lists 
the cell aspect ratio at the boundary, Ax/Ay (= b^) , for each case. The next 
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three columns show the number of iterations required to decrease the starting 
residual by 10~^ for three schemes previously discussed. Divergence is 
indicated by an entry "D." It is seen that scheme III is significantly less 
sensitive to grid stretching in the normal direction than are the 
discontinuous-a schemes, I and II, 


Table 1. 

Number of Iterations 

to Reduce 

Residual by 10 


Incompressible Circle 

Flow, 101 

by 21 Mesh 

Ax 

Ay 


Scheme 


I 

II 

III 

0.5 

44 

43 

34 

1 

72 

36 

51 

10 

68 

43 

47 

20 

99 

53 

48 

100 

212 

D 

34 

1000 

400 

D 

127 


As previously mentioned, the empirically-developed default settings in 
the TAIR code provide for an increased value of a near the surface; the 
default value satisfies the restriction (28) only for the first case in Table 
1, Ax/Ay = 0.5. For that case, convergence is obtained; the scheme diverges 
for the other listed cases for which the default setting violates restriction 
(28). In scheme I, the value of a near the surface met the restriction, and 
convergence was obtained for all the listed cases. 
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It should be noted again that the stability analysis presented herein is 

valid only for subsonic flow, when the AF2Y scheme is guaranteed to be 

hyperbolic in time. When the flow becomes locally supersonic, the linearized 

Eq. (1) will have a < 0, and a term which simulates d> must be added for 

stability [16], The effect of including such a term (e.g., in the second 

factor of Eq. (9)) has not been studied at present. With that cautionary 

remark, we present results for transonic cases in the next table. 

Table 2 presents results for two transonic cases for an NACA 0012 

airfoil: (1) Zero incidence with free-stream Mach number M = 0.85 and 

00 

(2) 2° incidence with M = 0.75. All cases were run with w - 1.8, but with 
different a- cycles. It can be seen that there is little difference in the 
convergence rate among the schemes, except that scheme II is noticeably slower 
than schemes I or III for case (2). 


Table 2. Number of Iterations to Decrease Residual by 10 ^ for 

Transonic Flow. NACA 0012, Default TAIR Mesh, 149 by 30 


Flow Condition 


Scheme 



I 

II 

III 

M = 0.85 

OO 

Zero incidence 

190 

174 

187 

M = 0.75 

OO 

2° incidence 

190 

360 

226 
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Three-Dimensional Version of AF2Y 

A three-dimensional (3-D) version of the AF2Y scheme is presented in Ref. 
7. It is different from the 2-D version discussed up to now, in that the 
factors are reversed in order. That is, the scheme can be written in the 
present context as: 


fa - ■£— 5 ) fb„ - — 5 ) (a - b. t ) A<{> .. = aaiLif)?. . + ab.fa - b.$ ) Ac{> , . , (35) 

^ b„ zz JK 2 a xx' 1 '- 1 y' Y jk£ T jk£ 2^ 1 y J r j,k-l,£ 


where 


L *jk£ = ^ a6 xx + 


b6 + c6 1 <j> . 
yy zz' T jk£ 


(36) 


Because the factors are reversed, we will refer to this scheme as AF2YR. 

Here the third coordinate direction is z, which can be thought of as the 
spanwise coordinate for a wing. The x- and y-coordinates are still the 
streamwise and normal coordinates as in the 2-D problem. The boundary 
operator corresponding to Eq. (8) is: 


Vj.l.t ■ ( a5 XX + 2bJ y + c5 zJ *3,1,1. " 2b4yS<]0 ' 


The scheme is implemented in three stages , as follows : 


1. (a - t— 6 ) g. = aa)L<j) a + ab~f . , . 

v b 2 zz' b j£ Y jk£ 2 j,k-l,J 


2 . 


( b 2 a ^xx^ f jk£ S j£ 


(37) 


(38) 

(39) 


3 


fa - b. $ ) Aft.,. = f.,„ 
1 y J Y jk£ jk£ 


(40) 
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The solution for f proceeds in planes, outward from the wing surface, using 
the tridiagonal Eqs. (38) and (39). The third stage (40) proceeds inward, 
solving for the correction in a bidiagonal march. 

Again, the main problem is how to initiate the first stage. In Reference 
7, the boundary condition used for f is 


f . 
J 


, 0 ,£ 


= 0. 


(41) 


We can again consider the more general boundary conditions studied previously, 


f j ,0 ,£ ~ 


(42) 


or 


f j,0,£ = bj ¥j,l ,£ 


(43) 


corresponding to Eqs. (18) and (30), respectively. Actually, condition (42) 
can only be approximately modeled in the 3-D problem, with some splitting 
error in the first two factors. That is, we can approximate Eq. (42) by 
solving, at k = 1: 


1. 


'BT “zz ; e j£ 


) g4» • 


B ¥ j,l,£ 


(44) 


2 . 


[ (1-Y) b„ - - 6 ]f . 
1 " n XX J •> 


L j , i , £ “ s ii' 


(45) 


Equation (43) is easily implemented by replacing a> in Eq. (38) by w + 8 
and setting f. n =0. 

J >*-) , x. 
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Stabllity Analysis of the 3-D AF2YR Scheme 

A VN analysis of the 3-D interior scheme shows that VN stability is 
achieved under restrictions (16) and (17). VN analysis of the boundary scheme 

(42) shows that sufficient conditions for stability of the VN boundary scheme 
are: 

0 < to < 1 - y (46) 

and 

Y < 1. (47) 

The same criteria are obtained in the 2-D counterpart of the AF2YR scheme with 
boundary condition (18). The corresponding criteria for boundary condition 

(43) are: 

0 < to + g < 1. (48) 

At this time we have no numerical experiments to test the stability and 
convergence of the 3-D boundary conditions (42) or (43) and the criteria (46) 
or (48). However, some comments about the use of AF2YR versus AF2Y are in 
order. 

In the AF2YR scheme, the use of boundary condition (42) or (43) makes the 
scheme parabolic at the surface; i.e., the time-like equation at the boundary 
is: 

o$ t = V 2 4> (49) 


where 


o = b 2 (1 -y)/u> 


(50) 
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for Eq. (42), and where 

a = b 2 /(w + 8) (51) 

for Eq. (43). In the case of AF2Y, the boundary equation remains hyperbolic , 
like the interior scheme, with 

" <f* yt = V 2 <j» (52) 

where 

a = b 2 /n. (53) 

It is felt that for this reason AF2Y may lead to faster convergence. It would 
appear that there is no difficulty in implementing such a scheme in 3-D, as: 


(“ + b lV V ' 

2. (ab 2 - cS J g jt - f ju + 



(54) 

(55) 

(56) 


The factors in the second and third stages could also be interchanged. The 
first stage is initiated by using Eq. (20), and the same stability and 
restrictions (26) and (27) hold. 
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CONCLUDING REMARKS 


We have studied the stability of the AF2Y scheme with several boundary 
conditions for the intermediate variable. The von Neumann method provides a 
useful tool for this study in view of the Goldberg-Tadmor theorem, and the 
results were verified in the two-dimensional case by the more complete GKSO 
theory. 

In general, the boundary schemes place a limitation on a and ui which 
is more restrictive than the requirements for the interior scheme. Since 
small a is desirable to damp low-frequency errors, one strategy involves 
increasing a at or near the boundary to meet the boundary restriction while 
using smaller o in the interior mesh. Such "discontinuous- a" schemes 
require further analysis of the stability at the line next to the discon- 
tinuity since the scheme there is no longer the interior scheme. They diverge 
on certain stretched grids. A safer strategy is to decrease w at the 
boundary to conform to the restrictions. This results in a more robust 
scheme; and it does not appear to suffer much, if any, loss in convergence 
rate. 

In regard to the 3-D AF2Y scheme, the current implementation in the TWING 
code involves a reversal of the factors from the 2-D TAIR code. We refer to 
this scheme as AF2YR. Although the reversal of the factors makes 
no difference in the interior (for the linear constant-coefficient case), 
there is a significant difference at the boundary. The AF2YR boundary scheme 
is parabolic in time as opposed to hyperbolic for AF2Y. For this reason, 
there may be a preference for the AF2Y, as in the TAIR code. 
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APPENDIX A 

Application of the GKSO Theory to the AF2Y Scheme 

In the GKSO theory [10], [11], the interior and boundary schemes are 
Considered as a coupled problem. Instead of substituting the Fourier 

solutions as in Eq. (22), the class of trial solutions is extended to 

,n _n i( jpAx) k , . , v 
<j)j k “ G e y (Al) 

where y is a complex number not restricted to lie on the unit circle in the 

complex plane. Fourier modes are retained in the direction tangential to the 

boundary under study. The trial solutions are substituted into the interior 
and boundary schemes, Eqs. (9) and (20), to obtain, respectively: 

[a + bj(l - i)] [- ctb 2 (y - 1) + 2A](G - 1) = ow[-2A + b(y - 2 + £■)] (A2) 

and 

[- ub 2 (y - 1) + 2A] (G - 1) = G[-2A - 2b(y - 1)], (A3) 

where Eq. (8) for . is used for the right-hand side of Eq. (A3) and 

B J .1 

where we have used the notation of Eq. (24). 

Equations (A2) and (A3) are two simultaneous equations for the unknowns 
G and y. In the theory, we are concerned only with values of y inside 
the unit circle, i.e», only those solutions which decay away from the 
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boundary. If the solution of Eqs. (A2) and (A3) yield G >1 for y < 1, 
the scheme is unstable. 

If Eq. (A3) is divided into Eq. (A2), G is eliminated; and there results 
an equation for y : 

ft[ay + bj(y - 1)][-2A + 2b (y - 1)] = aw[-2yA + b(y - 1)^]. (A4) 

First, it will be shown that for 

A = a( 1 - cos C) = 0 , 

there is a value of y inside the unit circle. When A = 0, Eq. (A4) reduces 
to two linear factors: 

(y “ l){[2ft (a + bj) - aa)]p - 2b^ ft + cuo} = 0. (A5) 

The root y = 1 is a solution of Eqs. (A2) and (A3) only when ft = m, 

corresponding to y = 1. (See Eq. (21).) Then 

G = 1 - a) (A6) 


and the restriction (16) must be satisfied. The other root is: 

2b^ft - am 

** 2ft (a + bj) - am 


(A7) 


which is less than 1.0 and is arbitrarily close to 1.0 as a approaches zero. 
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Using Eq. (A3), we can show that for any complex y such that its real 

part is less than 1, G <1 if and only if restrictions (26) and (27) are 

satisfied. Thus, let 

V = W R + iUj (A8) 

where y R and y^. are the real and imaginary parts. Substitution of Eq. 
(A8) into (A3) and multiplication by bj gives: 

[ab (1 - y R ) + 2Ab 1 - iabp I ](G - 1) = - 2nbjA + b (1 - y R ) - ibyj. (A9) 

2 

The condition G <1 then yields: 

n{A 2 b x (2 - n) + Ab(l - y R ) [a - + (2 - n)bj 

+ b 2 (a - fib^d - y R ) 2 + y 2 ]} > 0. (A10) 

For y < 1, the restrictions (26) and (27) are sufficient to ensure the 

inequality (A10) for arbitrary positive values of A, b^ , and b, regardless of 
the magnitude of y^. When A = 0, a value y <1 always occurs, as shown 

by Eq. (A7); and the scheme will be unstable unless restriction (27) is 
satisfied. Thus, restriction (27) is necessary; and when it is satisfied (for 
small a), restriction (26) will also be satisfied. 
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